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Abstract 

Ah initio determination of model Hamiltonian parameters for strongly correlated materials is a 
key issue in applying many-particle theoretical tools to real narrow-band materials. We propose a 
self-contained calculation scheme to construct, with an ah initio approach, and solve such a Hamil- 
tonian. The scheme uses a Wannier- function-basis set, with the Coulomb interaction parameter U 
obtained specifically for these Wannier functions via constrained Density functional theory (DFT) 
calculations. The Hamiltonian is solved by Dynamical Mean-Field Theory (DMFT) with the effec- 
tive impurity problem treated by the Quantum Monte Carlo (QMC) method. Our scheme is based 
on the pseudopotential plane-wave method, which makes it suitable for developments addressing 
the challenging problem of crystal structural relaxations and transformations due to correlation 
effects. We have applied our scheme to the "charge transfer insulator" material nickel oxide and 
demonstrate a good agreement with the experimental photoemission spectra. 
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I. INTRODUCTION 



The problem of inter-electron Coulomb interaction in solids is a many-particle problem 
which cannot be solved without major approximations. The most successful and widely 
used approach is electronic density functional theory (DFT)ii^ — within the local density 
approximation (LDA) or generalized gradient approximation (GGA), where all electrons feel 
the same one-particle potential defined by the electronic density distribution in the system. 
For wide-band materials, where the kinetic energy term dominates and the inter-electron 
Coulomb interaction can be treated in an average way with an energy (time) independent 
potential, this approach works very well. However for narrow-band materials this is not 
the case. The Coulomb interaction dominates over the kinetic energy, leading to strong 
correlations between electrons. Hence a static one-electron potential, as in DFT, becomes a 
bad approximation. Nevertheless, DFT can still give good results in such systems for some 
static integral properties. Moreover, DFT calculations can serve as a starting point for more 
sophisticated approaches designed to treat strongly correlated systems. 

The physics of strongly correlated systems was historically studied via solution of model 
Hamiltonians, such as the Hubbard^ and Anderson models^. In such Hamiltonians, the 
Coulomb interaction term is defined using a set of localized atomic-like orbitals centered on 
atomic sites. While the kinetic energy is invariant with the choice of the wave-functions- 
basis set of the model, the Coulomb interaction term significantly depends on the specific 
form of the atomic-like orbitals. 

Wannier functions^ are defined as Fourier transforms of Bloch functions, from wavevector 
to real space. They are considered nowaday as an optimal choice of basis set to construct 
model Hamiltonians,-!^ because they have a form of atomic centered localized orbitals and 
represent a complete basis set for the Bloch functions Hilbert space. However, Wannier 
functions are not uniquely defined and one needs to impose some additional conditions to 
make them unique. 

Another source of uncertainty in the construction of the Coulomb interaction term is 
the value of the onsite Coulomb repulsion parameter U. Attempts to simply determine 
it via the integral of the Coulomb potential multiplied by the squares of the orbital wave 
functions^'^ gave values 2-3 times larger than the experimental estimates. This is due to the 
neglect of strong relaxation and screening effects. The latter effects can be calculated using 
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perturbation theory, but the results strongly depend on the choice of the screening channels 
and the number of higher and lower lying states included in the expansion series. 

DFT can provide not only good data for the kinetic energy terms in the model Hamilto- 
nians, but also gives a practical alternative way to calculate the Coulomb parameter U using 
constrained DFT calculationsii'i^ii^'i^. In the constrained calculations, the DFT equations 
are solved with a fixed occupancy of the localized orbitals, and U is defined as a derivative 
of the orbital energy with respect to its occupancy. 

After its construction, the Hamiltonian for strongly correlated electrons needs to be 
solved in a way that should be as close to exact as possible. Recently Dynamical Mean- 
Field Theory (DMFT)^ii^ii^ii^ became a very popular tool to describe strongly correlated 
materials, especially when the numerically exact Quantum Monte Carlo (QMC)^^ method is 
used as a solver for the effective impurity problem. While neglecting, in its standard version, 
the inter-site correlations, DMFT fully accounts for the local dynamics. When the system 
is not very close to an ordering-disordering transition, single-site DMFT usually provides a 
satisfactory agreement between calculated and experimental spectra. 

Various methods are available to carry out DFT electronic structure calculations, which 
are based on different approximations for the wave-functions expansion. Any one of those 
methods may be used, in principle, to construct and solve the Hamiltonian for strongly 
correlated materials, if only the electronic spectral functions are needed as the results of the 
calculations. However, Coulomb correlation effects can lead to strong renormalization of the 
electron-lattice coupling and hence to complicated phase transitions due to lattice distor- 
tions. To describe such effects, lattice relaxation should be taken into account explicitly. 

The pseudopotential plane-wave method for first-principles DFT calculations is well 
suited and widely used to determine lattice relaxations in solids^i - in recent years this has 
also been extended to DFT plus onsite Coulomb interaction (DFT + U) calculations.— 
It would be desirable to use this method also as a basis to develop a self-contained scheme 
to determine, via DMFT and Wannier functions, the properties of strongly correlated ma- 
terials. So far, DMFT computations with Wannier functions have been implemented either 
using Muffin-Tin-Orbital (LMTO) methods, with Wannier functions constructed using the 
N-order muffin-tin-orbital scheme,— or more recently using some mixed-basis methods, with 
maximally localized Wannier functions,^ and some an adjustable U parameter value. 

In the present work, we propose a calculation scheme where all the terms of the Hamilto- 



3 



nian are generated, in a consistent way, within the ab initio pseudopotential plane-wave 
framework. Starting from the DFT pseudopotential-plane-waves method, we construct 
atomic-centered Wannier functions and produce the Hamiltonian kinetic-energy term in the 
basis of the Wannier functions. Then the value of the Coulomb parameter U for electrons in 
these Wannier functions is calculated via constrained DFT calculations. At the last stage, 
the DMFT-QMC method is used to solve the Hamiltonian. Recently, the electronic struc- 
ture of nickel oxide was successfully described within the DMFT method^^. In the current 
article NiO plays the role of a well known and already well described test system. The goal 
is to demonstrate the ability of our new method to reproduce the electronic structure of real 
system with a level of agreement comparable to previous LMTO-based DMFT calculations. 
Calculated spectral functions are compared with recent NiO DMFT-calculation results^ 
and with the experimental photoemission spectra. 

II. METHOD 

A. Wannier functions 

Wannier functions (WFs) \W!^) are defined as Fourier transforms of Bloch functions 

I^O^;^i:--^^l*->. (1) 

where T is the lattice translation vector, n the band number and k the reciprocal lattice 
vector. WFs are not uniquely defined because, in the single band case, there is a freedom of 
choice of the phases of the Bloch functions, I^E^nk)? as a function of k, and in the multiband 
case, any set of orthogonal linear combination of Bloch functions |\E'„k) could be used in 
([1]). The uncertainty in the WF's definition corresponds to a freedom of choice for a unitary 
transformation matrix U^^^ 

j 

One of the most commonly used approach to generate WFs was proposed by N. Marzari and 
D. Vanderbilt^. They use a condition of maximum localization of the WFs, that results in 
a variational procedure for the matrix Uj^^ . As an initial step before the variational process, 
a set of trial localized orbitals in the form of atomic orbitals was chosen and projected onto 
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the subspace of Bloch functions. Later-^ it was shown that this initial guess for the WFs 
of transition-metal oxides is usually so good that the variational procedure can be dropped 
and the projection of the trial orbitals onto the subspace of Bloch functions can be used to 
define the unitary transformation matrix . 

In the present work, we employ the pseudopotential method and a plane wave basis set. 
Hence, site centered pseudoatomic orbitals 0„ were chosen as a set of trial orbitals. 

Nonorthogonalized approximations to the WFs in the direct \W^) and reciprocal space 
\Wnk) are calculated as projection of the pseudoatomic orbitals onto a subspace of Bloch 
functions that is defined by setting an energy interval Ei < 5i(k) < E2 or some band 
numbers Ni < i < N2: 

\K) = E l^"k)e-^^^, (3) 

k 

\W„^)=J^\-^iu){-^i^.\M= Yl l^*k)(^.k|0„k). (4) 

i=Ni Ei<ei(k)<E2 

In the plane waves basis, Bloch functions and Bloch sums of pseudoatomic orbitals |0nk) = 
J2t ^^^'^\^n) '^^^ be decomposed as: 

|^^k) = 5^Q,q(k)|k-q), (5) 
q 

|0nk) = Ea„,q'(k)|k- q), (6) 

q' 

where n is a combined index jlma {j is the index for the atom, Im are the orbital and 
magnetic quantum numbers, respectively, and a is the spin projection). 

Using Eq. (j4]), together with these decompositions, one obtains: 

_ N2 N2 

IW^nk) = E l^.k)(^.k|0nk)= E |^.k)E<q(k)«.q'(k)(k-q|k-q') = 
i=Ni i=Ni q,q' 

N2 N2 ^ 

= E l*.k)Ec-q'(k)a,q'(k)= E &.,n(k)|*,k) = 

i=Ni q' i=Ni 

= E'^n,q(k)|k-q), (7) 

q 

N2 

bi,n(M) =Ycl^'(Myi,ci'{k), ^n,q(M) = (k)Ci,q(k) . (8) 

q' j=A''l 

To generate orthogonalized WFs, one should define the overlap matrix: 

0„„,(k) = (l^„k|W„'k) = Y Kni^)kn'i^)- (9) 

i=Ni 



Orthogonalized Wannier functions are then obtained as: 

\W:^) = Y^\Wn),)e^^^, where (10) 

k 

_ N2 

IW^nk) = $^(0„„.(k))-^|W^n'k) = ^^--(k)l*ik) = 5^C^„,q(k)|k- q) 

n' i=Ni q 

6,,„(k) = 5^(0„,,(k))-^6.^,(k) a;„,q(k) = 5^(0„„,(k))-^cu„,q(k). (11) 

n' n' 

For practical use — which includes generating the hamiltonian matrix for the kinetic energy 
term in the model Hamiltonian and performing constrained LDA/GGA calculation, one 
needs to compute the matrix elements of a given operator in the basis of the WFs. This can 
be conveniently done in reciprocal space. 

The matrix elements of the one-electron Hamiltonian in reciprocal space are defined as: 

/ N2 \ 

i7l^(k) = E l^.k)^.(k)(vi>,kl \Wm^) = 

\i=Ni J 
N2 

= E &r,„(k)6,^(k)£,(k), (12) 

i=Ni 

where ej(k) is the eigenvalue of the one-electron Hamiltonian for band i. 
In real space, the Hamiltonian matrix reads: 

hZI{t'-t) = {w:^\ (z E l^.k)e.(k)(vi/,k|) \w^') = 

\ k i=Ni J 

= E^nY(k)e^'^(^'-^). (13) 

k 

The Wannier functions occupancy matrix is given by: 

/ N2 \ 

\ k i=Ni / 

N2 

= E E blMkmim^^ik) - Ej), (14) 
k i=Ni 

where 6 is the step function and Ef is the Fermi energy. 

The transformation from plane waves to the Wannier functions basis is determined by 
Eqs. (ITl fTT]) and for the matrix elements by Eqs. f[T^ and f[T^ . It can be convenient also 
to back transform from WFs to the plane-waves basis. For example, if some constrained 
potential has a diagonal form in the WFs basis, Hnn' = SVnSnn', it can be written in the 
plane-waves basis as: 
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n 

H^onstr^^^ = (k-q|^,„„,,,(k)|k-q') = 

= J^(k-q|Vr„k)5K(W^nk|k-q') = 

n 

= ^<,q(k)5K^„,q,(k). (15) 

n 

B. Hamiltonian 

The localized Wannier functions are used as a basis set to define the Hamiltonian for 
strongly correlated materials: 



nm,TT' nm(nj^m),T 

Here c]jx(^"t) is the creation (annihilation) operator for an electron in the state defined by 
the Wannier function ri^x is the occupancy operator for this state, and -ffDC is a 

correction term to avoid double-counting for the Coulomb interaction that is already taken 
into account in DFT in an averaged way. 

The first Hamiltonian term, on the right-hand side of Eq. f|T6l) . corresponds to the kinetic 
energy, and the matrix elements H^^(T' — T) are defined by Eq. (JT3l) . The Coulomb in- 
teraction is described by the second term on the right-hand side of Eq. flTB]) . The Coulomb 
matrix Unm can be nonzero, in the most general case, for all orbitals. However, normally 
one takes into account only interactions between electrons in WFs having the symmetry of 
the transition-metal-ion d ot f orbitals. The corresponding matrix elements can be obtained 
using only two parameters: the direct Coulomb interaction parameter U and the exchange 



Coulomb interaction (Hund) parameter J (see for details Ref. l25l ). The double-counting 
correction is taken as: 

Hoc = ^eDcnnT- (17) 

riT 

Here we focus on a d-electron system and use for the double-counting potential: 

eDC = U{n,-^), (18) 



with Ud the number of d electrons^S. Eqs. ( ITTlfTSi) correspond to the assumption that the 
contribution to the DFT total energy from the Coulomb interaction between d electrons is 
given by: 

Edft = \und{nd - 1) (19) 

and the fact that one-electron energies in DFT are derivatives of the total energy with 
respect to the orbitals occupancy 

dE , , 

C. Constrained DFT calculations for the Coulomb interaction parameter U 

In a rigorous way, the Coulomb interaction parameter V for electrons in a state described 
by Wannier function Wnij) = {r\W^) should be calculated via the integral: 

U = J d\d'r'\Wn{T)\^U{r,r')\Wn{r')\^ (21) 

where the screened Coulomb interaction f/(r, r') is defined via the operator equation: 

U=[l-vP]-\ (22) 

with the bare Coulomb interaction, v{r, r') = l/(r — r'), and the polarization operator P: 

occ unocc f ^ _ -X 

P(r. r', = E E *(r)*-.-(r')*;(r)*,(r') [^^^^ - j;—^] ■ (^3) 



However various attempts^*^ to use Eqs. ( l2T]|23l) to calculate the Coulomb interaction pa- 
rameter U gave large variations in the resulting values. This is due to the many possibilities 
in choosing the channels of the screening via the definition of a set of occupied and unoccu- 
pied states in Eq. (l23l) . It either gives a strong underestimation of the U values, when one 
includes into the screening channels transitions between ci-states (that is clearly unphysical 
because electrons cannot screen themselves), or too large values of [/, when only states other 
than (i-states are included in the sums in Eq.f l23l) . but the summation over higher and lower 
energy states is not extended far enough. 

Another way to determine the Coulomb interaction parameter U is to use constrained 
DFT calculations,— 'i^iii^'^ where the screening and relaxation effects are taken into account 
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explicitly in a self-consistent procedure. If one uses the assumption that the contribution 
to the DFT total energy from the interaction between d-electrons is given by Eq. ffTQl) . then 
U can be calculated as a second derivative of the DFT total energy with respect to the 
(i-orbital occupancy: 

^~ d'n, ' ^^^^ 
using Eq. (l20!) this can be expressed via the first derivative of the one-electron energy e^: 

U = p^. (25) 

To use Eg. (1251) . one needs to perform DFT calculations with a constraint fixing the occu- 
pancy of (i-orbitals to certain values. In practice, this is done via an auxiliary potential in 
the form of a projection operator acting on (i-symmetry WFs IWJ): 

Hconstr = Y.\^n)SVn{W^\ (26) 
n 

(in reciprocal space this equation will take the form of Eq. lfTSll ). One-electron energy can be 
calculated as diagonal elements of the Hamiltonian matrix, in Eq. flT^ . and the corresponding 
occupancy as diagonal elements of the occupation matrix, in Eq. flT^ : 



= HZ^iO) (27) 

'I'd — ^dd 1 

and the derivative |^ (Eq. fl25l) ) is computed then numerically. 
D. Dynamical Mean-Field Theory 

The problem defined by the Hamiltonian (fT6l) can be solved by any of the methods 
developed to treat many-body effects. In the present work we have used Dynamical Mean- 
Field Theory (DMFT)^^'^^'^'^ which was recently found to be a powerful tool to numerically 
solve multiband Hubbard models. 

In DMFT, the lattice problem becomes an effective single-site problem, which has to be 
solved self-consistently for the matrix self-energy S and the local matrix Green function in 
the WFs basis set: 



Gnn'{e) =T-3— / dk 

VBZ 



{e-^i)l-H\\.)-ne)V\ , (28) 



where fi is the chemical potential, is the noninteracting part of the Hamiltonian f|T6|) : 

H' = J2 HnJiT' - T)c1tC„t' - Hoc (29) 

nm,TT' 

(in reciprocal space, the matrix elements for if''(k) can be calculated using Eq.( fT2l) ) and 
^{s) is the self-energy in the Wannier functions basis: 

^{e) =J2\Wn)^nn'{e){Wn'\. (30) 
nn' 

The DMFT single-site problem may be viewed as a self-consistent single-impurity Ander- 
son model^^. The corresponding local one-particle matrix Green function G can be written as 
a functional integral^ involving an action where the Hamiltonian of the correlation problem 
under investigation, including the interaction term with the Hubbard interaction, entera^^. 
The action depends on the bath matrix Green function Q through 

{g)-' = {Gr' + E. (31) 

To solve the functional integral of the effective single-impurity Anderson problem, various 
methods can be used: quantum Monte Carlo (QMC), numerical renormalization group 
(NRG), exact diagonalization (ED), noncrossing approximation (NCA), etc. (for a brief 



overview of the methods see Ref. |27l ). In the present work, the QMG^ method is used 
to solve the impurity problem. The real-frequencies single-particle spectral functions are 
computed using the maximum entropy method.— 



III. CALCULATIONS FOR NIO 

The calculation scheme described in Section Ull was applied in the present work to the 
problem of the classical Mott insulator nickel oxide, NiO. A DMFT study of the NiO elec- 
tronic structure was done recently^'' with LMTO-based Wannier functions and here we use 
previous work as a reference. Actually, in the Zaanen-Sawatzky-Allen classification,^- NiO 
is normally considered as a 'charge transfer insulator' where the minimal energy excitation 
across the gap happens between occupied oxygen p-states and unoccupied transition metal 
(i-states. This fact makes it crucial to include in the calculation not only partially filled 
(i-states (as is usually done for strongly correlated materials), but also explicitly take into 
account oxygen p-orbitals. 
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The pseudopotential plane- waves method, as implemented in the Quantum-ESPRESSO^- 
package, was used in the present work. The calculations were performed within the local den- 
sity approximation to density-functional theory, using Vanderbilt ultrasoft pseudopotential^ 
in the Rappe-Rabe-Kaxiras-Joannopoulos form.^° For the lattice parameter of NiO (rock- 
salt structure), we used the experimental value a = 4.193 A. A kinetic-energy cutoff of 40 
Ry was employed for the plane-wave expansion of the electronic states. The integrations in 
reciprocal space were performed using a (4,4,4) Monkhorst-Pack— k-point grid. 

The Wannier functions, introduced in Section [Tll are defined by the choice of the Bloch 
functions Hilbert space and by a set of trial localized orbitals. To show how this choice 
influences the results, we performed two calculations for the WFs. One with only bands 
formed by the Ni 3d states included in the Hilbert functional space and the second one 
where this space was extended by inclusion also of oxygen p-bands. In both cases, the 
projection was done on pseudoatomic wave functions. 

The results obtained without and with the p-bands are presented in Fig{T] and Fig|2l 
respectively. In both cases, the ci-bands obtained in the pseudopotential calculations are 
exactly reproduced by the bands obtained diagonalizing the k-space Hamiltonian matrix 
(Eq. f|T2l) ). However the spatial distributions of the corresponding Wannier functions are 
different in the two cases. To demonstrate that, spatial distribution of (ia^^-like Wannier 
function charge density was plotted (see Fig {1] and Fig|2]). Omission of the oxygen bands in 
the first calculation results in the appearance of a contribution from p orbitals on neighboring 
oxygen atoms in the Wannier function obtained by projection of the Ni-atom 3d orbitals 
(see FigUI). Extending the Hilbert space with explicit inclusion of oxygen p bands results in 
a Wannier function that is nearly a pure rf-orbital (see Figj2]). In the following, we will use 
WF obtained with the full sets of bands and pseudo-atomic orbitals including all Ni 3d and 
O 2p states. 

The partial densities of states calculated from the WF hamiltonian are shown in Fig|3l 
One can distinguish the occupied oxygen p-band and, at higher energy, the partially occupied 
d-hand. The t2g sub-band is nearly filled, while the Fermi level is located essentially in the 
middle of the Cg sub-band. 

The calculation of Coulomb interaction parameter U was done using straightforward 
constrained DFT procedure^. We performed a self-consistent calculation in a supercell 
containing two formula units. The constrained potential 6Vn (Eq. fl26p ) on the first Ni atom 
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FIG. 1: (Color online) Top frame: Calculated LDA band structure of NiO (solid line) and bands 
obtained from the five d-type Wannier functions (circles). The energy window used to compute 
the Wannier functions was set to [-2. 5, +1.5] eV from the Fermi level. Bottom frame: Charge 
density of the corresponding d^z-like Wannier function (plotted using the XCrysDen^ package). 

was positive and equal to 0.1 Ry, and for the second Ni atom it had the same magnitude and 
opposite sign. The Wannier functions occupations and energies were calculated and the U 
value was then evaluated as described in Sec. [Tll using Eqs. (12511271) . The resulting value for 
the Coulomb interaction parameter U was found to be 6.6 eV. While this is smaller than the 
value of 8 eV obtained from constrained LDA calculations for linear-muffin-tin orbitals,— it 
is larger than the value of of 4.5 eV obtained from linear- response-theory calculations for 
pseudo-atomic orbitals.— 

The interacting Hamiltonian, constructed from the pseudopotential plane-waves calcula- 
tions, was solved by DMFT-QMC simulations. In the QMC simulations, the inverse tem- 
perature value was P = 10 eV~^ (T = 1160i^), and we used 80 time slices and 4 x 10^ 
Monte Carlo sweeps. Other DMFT+QMC-calculation details are the same as in^. The 
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FIG. 2: (Color online) Top frame: Calculated LDA band structure of NiO (solid line) and bands 
obtained from the eight Wannier functions (circles). The energy window used to compute the 
Wannier functions was set to [-8. 5, +1.5] eV from the Fermi level. Bottom frame: Charge density 
of the corresponding d^^-like Wannier function. 

DMFT calculations result in strong changes in the spectral functions (see FigHj) compared 
to the LDA densities of states (Fig|3]), consistent with recent DMFT calculations based on 
the LMTO method.—*^ First of all, the correlated spectra show a large-gap insulator in 
contrast to the metallic LDA solution. The half filled Cg band is now split into an empty 
upper-Hubbard band and an occupied lower-Hubbard band. In addition, with explicit in- 
clusion of the oxygen 2p band, the lower Hubbard Cg band is found to largely overlap in 
energy with the oxygen p-band and to strongly hybridize with it. Furthermore, although 
the t2g band remains occupied, it is shifted down in energy so that it also overlaps with the 
oxygen band. 

In FigHl we present a comparison of the experimental photoemission (direct and inverse) 
spectrar^ with the spectral functions obtained from our DMFT pseudopotential calculations. 
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FIG. 3: (Color online) Orbital resolved density of states for the non-interacting LDA hamiltonian. 
Solid black line: Cg states of Ni, dashed red line: t2g states of Ni, solid green line: p states of O. 
Here and later Fermi energy corresponds to zero. 

The photoemission spectra measured with 120 eV photons reflect more Sc? Ni states and those 
with 66 eV have more oxygen contributions, due to the difference in the cross section values. 
One can observe a good agreement between the experimental and calculated spectra. The 
calculations reproduce a large band gap as well as remarkable redistribution of spectral 
weight from the top of the valence band to lower energies, going from the 120 eV to the 66 
eV spectrum. Comparison with the theoretical orbital-resolved result in FigHJ indicate that 
this is consistent with a decreasing (i-states spectral density contribution and an increasing 
p-states contribution to the spectrum. 

IV. CONCLUSION 

We have presented a computational scheme for strongly correlated materials where, in the 
framework of the pseudopotential plane-waves method, atomic-centered Wannier functions 
are calculated and the Hamiltonian matrix elements are evaluated in the Wannier-functions 
basis. Then, via constrain DFT calculations, the Coulomb interaction parameter U is ob- 
tained for electrons in these Wannier functions, allowing us to generate in a consistent way 
all necessary components of the Hamiltonian. This Hamiltonian is solved by Dynamical- 
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FIG. 4: (Color online) Orbital resolved density of states of NiO for the interacting Hamiltonian. 
Solid black line: Cg states of Ni, dashed red line: t2g states of Ni, solid green line: p-states of O. 

Mean-Field Theory with the numerically exact Quantum Monte-Carlo method as effective- 
impurity-problem solver. We have applied this scheme to nickel oxide and demonstrate good 
agreement with the experimental photoemission data. 
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